Take-home Exercise 1

My Data Processing Journey 🛠️

I started running the rendering process 3 days ago.
I encountered multiple memory issues, missing packages, and code debugging.
Every time I tried rendering, it took 5 hours to finish… or crash! 😅

After numerous debugging sessions, I managed to reduce the data size, clean up the code, and finally, after 3 days, the rendering is successful! 🎉

What I Learned

Even though I shrank the data to 1/10 of its original size and optimized the code and follow step in piazza, rendering was still a challenge.
But the process was enlightening, and I have gained valuable debugging skills along the way.

---

Current Status 🏃‍♂️

Good news: The code is now running smoothly.
I’m currently working on the K-function analysis, as the next milestone.

Thank you for your patience, and I’ll keep you updated on my progress! 😊

Setting the Scene

As urbanization increases globally, cities are becoming more congested with vehicles, leading to a higher incidence of traffic accidents. Effective urban planning and road safety interventions rely on a deep understanding of where and when accidents are most likely to occur. The ability to identify traffic hotspots and patterns over time is critical for city authorities and policymakers to develop more targeted safety measures.

In this study, we examine the traffic accident data for the Bangkok Metropolitan Region (BMR), one of Southeast Asia’s largest urban conglomerations. The analysis focuses on spatial and temporal dimensions to identify accident hotspots across different road networks. Using Kernel Density Estimation (KDE) and Network Constrained K-function analysis, we aim to detect areas of high risk and analyze how accident patterns vary across time periods (e.g., day vs. night, weekdays vs. weekends).

This project builds on advanced geospatial analysis techniques and is designed to provide critical insights into the spatial concentration of traffic accidents. The results will help decision-makers identify the most accident-prone areas and develop interventions to reduce the risk of accidents in these hotspots.

Objectives

The main objective of this take-home exercise is to explore spatio-temporal dynamics and identify traffic accident hotspots in the Bangkok Metropolitan Region (BMR). The specific goals are:

Accident Hotspot Detection: Use Kernel Density Estimation (KDE) to visualize and detect areas of high accident density across the BMR road network. This will be done both at a regional level and for selected high-risk provinces (e.g., Bangkok and Samut Prakan).

Spatio-Temporal Analysis: Perform temporal analysis by splitting the data into different time periods (e.g., day vs. night, weekdays vs. weekends) to understand how accident patterns evolve over time.

Network Constrained Analysis: Apply Network Constrained Kernel Density Estimation (NKDE) and K-function analysis to detect clustering patterns of accidents on road networks. This analysis will help understand how accidents are distributed along specific road networks and identify clusters of incidents.

Visualization and Reporting: Use geospatial visualization tools to effectively communicate findings. This includes generating interactive heatmaps and visual comparisons of accident density for different time periods and road networks.

The Data

Road Accident Data (2019-2022): This dataset includes information on over 12,986 road accidents that occurred in the BMR region. Each entry contains detailed information such as the time and location of the incident, the type of road, and other contextual factors. The key attributes in this dataset include:

incident_datetime: The date and time when the accident occurred. province: The province where the accident took place. road_type: Information about the road category where the accident occurred (e.g., highways, secondary roads). Bangkok Metropolitan Region Roads (from OpenStreetMap): The road network data includes detailed geographical information on the roads in the BMR. The dataset provides the structure of roads as LINESTRING geometries and includes information on road types, length, and geometry details.

Thailand Subnational Administrative Boundaries: This dataset includes the administrative boundaries of the provinces in the BMR. It allows for an understanding of accident distributions across different administrative regions (e.g., Bangkok, Samut Prakan, Pathum Thani, etc.).

The combination of these datasets enables us to perform detailed geospatial analysis, including network-based density estimations and spatial clustering detection. The accident data is pre-processed to clean any inconsistencies and transformed to appropriate coordinate reference systems (CRS) to align with road network data.

Getting Started

To begin, several R packages are loaded for different purposes such as spatial data processing, network analysis, and visualization:

Code

pacman::p_load(sf, sfdep, tmap, tidyverse, spNetwork, spacetime, spatstat, ggplot2)

sf (Simple Features): Provides tools for manipulating and analyzing spatial data, specifically through sf objects which store geometries like points, lines, and polygons.

sfdep: Used for spatial dependence analysis, including functions for working with neighbors and performing spatial lag models. tmap: A package designed for creating both static and interactive maps. It offers tools for geospatial data visualization, which is critical for this analysis to present accident hotspots and spatial patterns effectively. tidyverse: A collection of R packages that simplifies data manipulation, visualization, and exploration through functions like dplyr, ggplot2, and tibble.

spNetwork: A specialized package for network-based spatial analysis. It supports functions like network-constrained Kernel Density Estimation (NKDE) and network-constrained K-function analysis.

spacetime: Provides classes and methods for spatio-temporal data, allowing analysis over both space and time. stpp: Designed for space-time point process modeling and simulation. Useful for detecting spatio-temporal clusters and patterns in the accident data.

spatstat: A package for spatial point pattern analysis, providing tools for calculating K-functions and conducting spatial statistics over road networks.

##Importing Data ### 1. Load and Clean Accident Data

Code

car_acc <- read_csv("data/rawdata/thai_road_accident_2019_2022.csv") %>%

   # Remove rows with missing longitude or latitude
  filter(!is.na(longitude) & !is.na(latitude)) %>%

  # Create new columns for month and day of the week
  mutate(Month_num= month(incident_datetime)) %>%
  mutate(Month_fac= month(incident_datetime,label=TRUE,abbr=TRUE))%>%
  mutate(dayofweek=day(incident_datetime))%>%
  # Convert the data frame to a spatial sf object
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%

  # Transform to a different CRS (Coordinate Reference System)
  st_transform(crs = 32647)

read_csv(): Reads the CSV file containing the road accident data. filter(!is.na(longitude) & !is.na(latitude)): Removes any rows where either longitude or latitude is missing, ensuring valid spatial data. mutate(Month_num = month(incident_datetime)): Creates a new column Month_num that extracts the month from the incident_datetime column. mutate(Month_fac = month(incident_datetime, label = TRUE, abbr = TRUE)): Creates a new column Month_fac that formats the month as a labeled and abbreviated factor. mutate(dayofweek = day(incident_datetime)): Creates a new column dayofweek that extracts the day of the month from incident_datetime. st_as_sf(coords = c(“longitude”, “latitude”), crs = 4326): Converts the data into a spatial object with coordinates based on longitude and latitude. CRS 4326 is the WGS84 geographic coordinate system. st_transform(crs = 32647): Transforms the spatial object to UTM zone 47N (EPSG: 32647), which is often used for geographic data in Thailand.

2. Filter Accident Data for BMR (Bangkok Metropolitan Region)

Code

 bmr_acc <- car_acc %>%
   filter(province_en %in% c("Bangkok", "Nonthaburi", "Nakhon Pathom", "Pathum Thani", "Samut Prakan", "Samut Sakhon"))

filter(province_en %in% c(…)): Filters the accident data to include only accidents occurring in the Bangkok Metropolitan Region, i.e., the six listed provinces.

3. Save the Filtered Data for BMR

Code

write_rds(bmr_acc,"data/rds/bmr_acc.rds")

4. Load Administrative Boundaries

Code

boundaries <- st_read(dsn = "data/rawdata",
                        layer = "tha_admbnda_adm1_rtsd_20220121")

5. Filter Boundaries for BMR

Code

bmr_boundary <- boundaries %>%
  filter(ADM1_EN %in% c("Bangkok", "Nonthaburi", "Nakhon Pathom", "Pathum Thani", "Samut Prakan", "Samut Sakhon"))

6. Save Filtered BMR Boundaries

Code

write_rds(bmr_boundary, "data/rds/bmr.rds")

7. Load Road Network Data

Code

roads <- st_read(dsn = "data/rawdata",
                        layer = "hotosm_tha_roads_lines_shp")

8. Assign CRS to Road Network

Code

st_crs(roads) <- 4326

9. Clip Roads to BMR Boundaries

Code

bmr_roads <- st_intersection(roads, bmr_boundary)

10. Save the Filtered BMR Roads Data

Code

saveRDS(bmr_roads,"data/rds/bmr_roads.rds")

11. Load the Date from RDS

Code

bmr <- read_rds("data/rds/bmr.rds")

Code

bmr_roads <- read_rds("data/rds/bmr_roads.rds")

Code

bmr_acc <- read_rds("data/rds/bmr_acc.rds")

EDA (Exploratory Data Analysis)

The EDA focuses on visualizing the basic geographic features and accident points within the Bangkok Metropolitan Region (BMR). Here’s a breakdown of each component of the code and the result:

EDA 1: Visualization of Traffic Accidents in the BMR

Code

tm_shape(bmr) + 
  tm_polygons() + 
  tm_shape(bmr_roads) + 
  tm_lines() +  
  tm_shape(bmr_acc) + 
  tm_dots(col = "red", size = 0.1) +  
  tm_layout(title = "Road Traffic Accidents in BMR")

tmap_mode("plot")

Explanation of the Code: tm_shape(bmr) + tm_polygons():

This part of the code visualizes the boundary of the BMR region by using the polygon data from bmr. The boundary defines the geographic limits of the region under analysis. tm_shape(bmr_roads) + tm_lines():

This adds the road network to the visualization. The roads are crucial for understanding the spatial relationship between the accident locations and the infrastructure in the region. tm_shape(bmr_acc) + tm_dots(col = “red”, size = 0.1):

This plots the accident locations on the map as red dots. The small dot size ensures that even when multiple accidents occur close to each other, the map remains readable. tm_layout(title = “Road Traffic Accidents in BMR”):

This switches the map into interactive mode, allowing users to zoom, pan, and explore the geographic data more dynamically. Result Analysis: Accident Distribution: The red dots represent individual traffic accidents in the BMR. The visualization clearly shows that the accidents are concentrated along the major roads and highways, suggesting that higher traffic volumes and more complex intersections may be contributing to the occurrence of accidents. High-risk Areas: Some road segments appear to have a dense cluster of accidents, which could indicate accident-prone areas. These could be major intersections, busy highways, or areas with challenging road conditions.

EDA 2: Monthly Distribution of Traffic Accidents in BMR

Code

monthly_accidents <- bmr_acc %>%
  group_by(Month_fac) %>%
  summarise(count = n())

tm_shape(bmr) +
  tm_borders() +  
  tm_shape(bmr_acc) +
  tm_dots(col = "Month_fac", palette = "Reds", size = 0.05) +  
  tm_facets(by = "Month_fac", free.coords = FALSE) + 
  tm_layout(title = "Accidents by Month in BMR")

Code

tmap_mode("plot") 

This visualization illustrates the monthly distribution of traffic accidents in the Bangkok Metropolitan Region (BMR). Accidents are grouped and color-coded by month, showing how accident occurrences vary across the calendar year. This helps in identifying temporal patterns and trends related to seasonal factors, road usage, or external influences (e.g., weather, holidays) that may contribute to the increase or decrease in accidents during specific months. The map uses a facet approach, allowing for a detailed month-by-month comparison, making it easier to visually assess changes in accident distribution across different time periods.

EDA 3: Heatmap of Traffic Accident Density in BMR Provinces

Code

bmr$ADM1_EN <- trimws(tolower(bmr$ADM1_EN))
bmr_acc$province_en <- trimws(tolower(bmr_acc$province_en))


setdiff(bmr$ADM1_EN, bmr_acc$province_en)

character(0)

Code

setdiff(bmr_acc$province_en, bmr$ADM1_EN)

character(0)

Code

accidents_by_province <- bmr_acc %>%
  group_by(province_en) %>%
  summarise(accident_count = n())

bmr_df <- st_drop_geometry(bmr)


accidents_by_province <- accidents_by_province %>% 
  mutate(province_en = tolower(province_en)) 
bmr_df <- bmr_df %>% 
  mutate(ADM1_EN = tolower(ADM1_EN))

bmr_with_accidents <- left_join(bmr_df, accidents_by_province, by = c("ADM1_EN" = "province_en"))
# Reorder columns to ensure ADM1_EN is the first
bmr_with_accidents <- bmr_with_accidents %>%
  select(ADM1_EN, everything())


bmr_with_accidents_sf <- st_as_sf(bmr_with_accidents, geometry = bmr$geometry)
tmap_mode("view")
tm_shape(bmr_with_accidents_sf) +
  tm_polygons("accident_count", 
              style = "pretty", 
              palette = "Reds", 
              title = "Accidents per Province", 
              popup.vars = c("Province" = "ADM1_EN", "Accidents" = "accident_count"), 
              popup.format = list(digits = 0)) + 
  tm_layout(title = "Accident Count Heatmap in BMR", legend.format = list(text.separator = "-"))

Explanation: Data Preparation:

Accident Data Processing: We started by reading in the accident dataset and filtering out rows where the coordinates (latitude, longitude) were missing. We then extracted the relevant date components (month and day of the week) from the accident timestamps for future analysis. The dataset was then converted into a spatial object using the appropriate geographic coordinate system (CRS 4326). Boundary and Road Data: Similarly, we read in the boundary and road data for the entire country of Thailand. We filtered out only the relevant provinces that belong to the Bangkok Metropolitan Region (BMR), which includes Bangkok, Nonthaburi, Nakhon Pathom, Pathum Thani, Samut Prakan, and Samut Sakhon. The roads were intersected with the BMR boundary to limit them to the region we are interested in. Province-Level Accident Count:

Standardization of Names: Before merging the datasets, we ensured consistency in the naming of provinces across the two datasets (bmr and bmr_acc). We converted all province names to lowercase and removed any extra spaces. Accident Count: We then grouped the accident data by province and counted the total number of accidents for each province. Merging Data:

We merged the accident counts with the BMR boundary data to assign the accident data to the spatial geometry of each province. This allows us to visually analyze accident density within each province. Visualization:

We used the tmap library to create an interactive heatmap. The accident density for each province is color-coded in shades of red, where darker shades indicate higher accident counts. Hovering over a province will display the exact count of accidents for that region. Result Analysis: The heatmap generated from the data analysis allows us to visually assess which provinces in the BMR have higher accident densities. The interactive aspect lets users hover over each province to see the accident count, making it easy to compare different regions. Key observations from this visualization could include:

Bangkok likely has the highest density of accidents due to its high population density and traffic volume. Samut Prakan and Nonthaburi may also show significant accident counts due to their proximity to Bangkok and the presence of major highways. Further analysis can be done by looking into different time periods (day vs. night, weekdays vs. weekends) or by exploring specific road networks using KDE methods to identify high-risk areas. The insights gained from this analysis can help in urban planning, improving road safety measures, and optimizing traffic management to reduce accidents in high-density areas.

Network KDE (NKDE) Analysis

1.road-data preparation：

Code

bmr_roads <- bmr_roads %>%
  filter(highway %in% c("motorway", "trunk", "primary", "secondary"))

bmr_city <- st_read(dsn = "data/rawdata",
                        layer = "tha_admbnda_adm2_rtsd_20220121")

Reading layer `tha_admbnda_adm2_rtsd_20220121' from data source 
  `/Users/lucasluo/Desktop/SMU/Courses/Term3 Aug-Dec/ISSS626-Geospatial Analytics and Applications/lucasluo6/ISSS626/Take-home_Ex/Take-home_Ex01/data/rawdata' 
  using driver `ESRI Shapefile'
Simple feature collection with 928 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 97.34336 ymin: 5.613038 xmax: 105.637 ymax: 20.46507
Geodetic CRS:  WGS 84

Code

# Step 1: 
SamutPrakan_roads <- bmr_roads %>%
  filter(ADM1_EN == "Samut Prakan")

bangkok_roads <- bmr_roads %>%
  filter(ADM1_EN == "Bangkok")

# Step 2: 
SamutPrakan_boundary <-bmr_city %>%
  filter(ADM1_EN == "Samut Prakan")

bangkok_boundary <- bmr_city %>%
  filter(ADM1_EN == "Bangkok")


SamutPrakan_roads_intersection <- st_intersection(SamutPrakan_roads,SamutPrakan_boundary)


bangkok_roads_intersection <- st_intersection(bangkok_roads, bangkok_boundary)

2.Choosing city-data：

Code

bangkok_boundary <- bmr_city %>%
  filter(ADM1_EN == "Bangkok")

SamutPrakan_boundary <- bmr_city %>%
  filter(ADM1_EN == "Samut Prakan")

3.Choosing acc-data：

Code

bangkok_acc_data <- bmr_acc %>%
  filter(province_en == "bangkok")

SamutPrakan_acc_data <- bmr_acc %>%
  filter(province_en == "samut prakan")

4.linestring：

Code

target_crs <- 32647  # UTM Zone 47N

bangkok_acc_data <- st_transform(bangkok_acc_data, crs = target_crs)
bangkok_roads_intersection <- st_transform(bangkok_roads_intersection, crs = target_crs)
bangkok_boundary <- st_transform(bangkok_boundary, crs = target_crs)


acc_in_bangkok <- st_intersection(bangkok_acc_data, bangkok_boundary)
roads_in_bangkok <- st_intersection(bangkok_roads_intersection, bangkok_boundary)


roads_in_bangkok <- roads_in_bangkok %>%
  filter(st_geometry_type(roads_in_bangkok) %in% c("LINESTRING", "MULTILINESTRING"))


roads_in_bangkok <- st_cast(roads_in_bangkok, "LINESTRING", group_or_split = TRUE)


SamutPrakan_acc_data <- st_transform(SamutPrakan_acc_data, crs = target_crs)
SamutPrakan_roads_intersection <- st_transform(SamutPrakan_roads_intersection, crs = target_crs)
SamutPrakan_boundary <- st_transform(SamutPrakan_boundary, crs = target_crs)

acc_in_SamutPrakan <- st_intersection(SamutPrakan_acc_data, SamutPrakan_boundary)
roads_in_SamutPrakan <- st_intersection(SamutPrakan_roads_intersection, SamutPrakan_boundary)


roads_in_SamutPrakan <- roads_in_SamutPrakan %>%
  filter(st_geometry_type(roads_in_SamutPrakan) %in% c("LINESTRING", "MULTILINESTRING"))

roads_in_SamutPrakan <- st_cast(roads_in_SamutPrakan, "LINESTRING", group_or_split = TRUE)

5.acc_in_city：

Code

target_crs <- 32647  

bmr_acc_data <- st_transform(bmr_acc, crs = target_crs)
bmr_boundary <- st_transform(bmr, crs = target_crs)

acc_in_bmr <- st_intersection(bmr_acc_data, bmr_boundary)

roads_in_bangkok_lines <- st_cast(roads_in_bangkok, "LINESTRING")
roads_in_SamutPrakan_lines <- st_cast(roads_in_SamutPrakan, "LINESTRING")

Code

lixels_bangkok <- lixelize_lines(roads_in_bangkok_lines,
                         10000,       
                         mindist = 5000)  


samples_bangkok <- lines_center(lixels_bangkok)

Code

roads_in_SamutPrakan_lines <- st_cast(SamutPrakan_roads_intersection, "LINESTRING")

lixels_SamutPrakan <- lixelize_lines(roads_in_SamutPrakan_lines,
                                      10000,        
                                      mindist = 5000)  

samples_SamutPrakan <- lines_center(lixels_SamutPrakan)

6.NKDE:

Code

acc_in_bangkok <- st_as_sf(bangkok_acc_data)

acc_in_SamutPrakan <- st_as_sf(SamutPrakan_acc_data)

nkde_result_bangkok <- nkde(
  lines = lixels_bangkok,                     
  events = acc_in_bangkok,                     
  w = rep(1, nrow(acc_in_bangkok)),            
  samples = samples_bangkok,                   
  kernel_name = "quartic",                     
  bw = 500,                                    
  div = "bw",                                 
  method = "simple",                          
  grid_shape = c(200, 200),                    
  verbose = FALSE                               
)

nkde_result_SamutPrakan <- nkde(
  lines = lixels_SamutPrakan,                 
  events = acc_in_SamutPrakan,                
  w = rep(1, nrow(acc_in_SamutPrakan)),      
  samples = samples_SamutPrakan,              
  kernel_name = "quartic",                     
  bw = 500,                                   
  div = "bw",                                 
  method = "simple",                          
  grid_shape = c(200, 200),                   
  verbose = FALSE                               
)

7.Density：

Code

samples_bangkok$density <- nkde_result_bangkok
lixels_bangkok$density <- nkde_result_bangkok

samples_bangkok$density <- samples_bangkok$density * 10000
lixels_bangkok$density <- lixels_bangkok$density * 10000

samples_SamutPrakan$density <- nkde_result_SamutPrakan
lixels_SamutPrakan$density <- nkde_result_SamutPrakan

samples_SamutPrakan$density <- samples_SamutPrakan$density * 10000
lixels_SamutPrakan$density <- lixels_SamutPrakan$density * 10000

8.Map：

Code

tmap_mode('view')
# Bangkok
tm_shape(lixels_bangkok) +
  tm_lines(col = "density", palette = "-RdYlBu", title.col = "Density", lwd = 5, breaks = c(0, 0.1, 0.2, 0.3, 0.4, 0.5)) + 
  tm_shape(acc_in_bangkok) +
  tm_dots(size = 0.1, col = "blue", alpha = 0.5, title = "Accidents")

Code

# SamutPrakan
tm_shape(lixels_SamutPrakan) +
  tm_lines(col = "density", palette = "-RdYlBu", title.col = "Density", lwd = 5, breaks = c(0, 0.1, 0.2, 0.3, 0.4, 0.5)) + 
  tm_shape(acc_in_SamutPrakan) +
  tm_dots(size = 0.1, col = "blue", alpha = 0.5, title = "Accidents")

##Perform K-function analysis

Spatio-temporal analysis

Step 1: Accident Temporal Patterns Aggregate Accident Data by Time: Extract temporal information from the accident data, such as hour of the day, day of the week, or month. This will help you see if accidents tend to happen during certain times of the day or specific months. Visualize Temporal Patterns: Use time series plots or heatmaps to visualize how accident frequencies change over time. You can also look for seasonality or trends (e.g., higher accident rates during specific months or on weekends).

Code

# Extract the hour, day, and month from the accident timestamps
bangkok_acc_data$hour <- lubridate::hour(bangkok_acc_data$incident_datetime)
bangkok_acc_data$day <- lubridate::wday(bangkok_acc_data$incident_datetime, label = TRUE)
bangkok_acc_data$month <- lubridate::month(bangkok_acc_data$incident_datetime, label = TRUE)

# Plot accident frequencies by hour of the day
ggplot(bangkok_acc_data, aes(x = hour)) +
  geom_histogram(binwidth = 1, fill = "steelblue", color = "black") +
  labs(title = "Accidents by Hour of the Day", x = "Hour", y = "Accident Count")

Code

# Plot accident frequencies by day of the week
ggplot(bangkok_acc_data, aes(x = day)) +
  geom_bar(fill = "darkorange", color = "black") +
  labs(title = "Accidents by Day of the Week", x = "Day", y = "Accident Count")

Code

# Plot accident frequencies by month
ggplot(bangkok_acc_data, aes(x = month)) +
  geom_bar(fill = "darkgreen", color = "black") +
  labs(title = "Accidents by Month", x = "Month", y = "Accident Count")

Step 2: Accident Hotspot Analysis Over Time Group Accidents by Time Intervals: split the accident data into different time intervals (e.g., hourly, daily, or monthly), and then analyze whether certain hours or months are more prone toaccidents. Overlay Spatial Data with Time: Overlay the accident data on the road network for different time periods to visualize whether spatial hotspots change at certain times of day, such as during rush hours.

Code

# Group accidents by time intervals (hourly, daily, etc.)
hourly_acc <- bangkok_acc_data %>%
  group_by(hour) %>%
  summarize(accident_count = n())

# Overlay time-based accident data with road data (for a specific time window)
rush_hour_acc <- bangkok_acc_data %>% filter(hour >= 7 & hour <= 9)  # Example: Rush hour 7-9 AM

# Use tmap to create spatial visualizations for rush hour accidents
tmap_mode('plot')
tm_shape(bangkok_roads) +
  tm_lines(col = "gray", lwd = 1) +
  tm_shape(rush_hour_acc) +
  tm_dots(size = 0.1, col = "red", alpha = 0.7) +
  tm_layout(title = "Accidents during Rush Hour (7-9 AM)")

Step 3: Spatio-Temporal Hotspot Detection Use Kernel Density Estimation (KDE) for Different Time Periods:

Calculate KDE for accident density, but now do it separately for different time periods (e.g., day vs. night, weekdays vs. weekends). This allows us to see how the spatial distribution of accidents changes over time. Network-based KDE by Time Intervals: ’’’’’’’

Code

# Set target CRS as UTM Zone 47N (EPSG: 32647)
target_crs <- 32647

# Transform the roads to UTM (projected coordinate system)
bk_roads_projected <- st_transform(bangkok_roads, crs = target_crs)

# Ensure simple geometries (convert MULTILINESTRING to LINESTRING if necessary)
bmr_roads_simple <- st_cast(bk_roads_projected, "LINESTRING")

# Transform accident data to UTM as well
rush_hour_acc_projected <- st_transform(rush_hour_acc, crs = target_crs)

# Convert time to hour (0-23)
bangkok_roads$hour <- as.integer(format(bangkok_roads$incident_datetime, "%H"))

# Convert date to day of the week (1 = Monday, 7 = Sunday)
bangkok_roads$weekday <- as.integer(format(bangkok_roads$incident_datetime, "%u"))

# Convert time to hour (0-23)
bangkok_acc_data$hour <- as.integer(format(bangkok_acc_data$incident_datetime, "%H"))

# Convert date to day of the week (1 = Monday, 7 = Sunday)
bangkok_acc_data$weekday <- as.integer(format(bangkok_acc_data$incident_datetime, "%u"))

# Define daytime and nighttime accidents (e.g., daytime: 6 AM to 6 PM)
daytime_acc <- bangkok_acc_data%>% filter(hour >= 6 & hour < 18)
nighttime_acc <- bangkok_acc_data %>% filter(hour < 6 | hour >= 18)

# Define weekday and weekend accidents (e.g., weekdays: Monday to Friday)
weekday_acc <- bangkok_acc_data %>% filter(weekday >= 1 & weekday <= 5)
weekend_acc <- bangkok_acc_data %>% filter(weekday == 6 | weekday == 7)

# Ensure your roads and accidents are in the projected CRS (UTM)
daytime_acc_projected <- st_transform(daytime_acc, crs = target_crs)
nighttime_acc_projected <- st_transform(nighttime_acc, crs = target_crs)
weekday_acc_projected <- st_transform(weekday_acc, crs = target_crs)
weekend_acc_projected <- st_transform(weekend_acc, crs = target_crs)

Code

# Perform NKDE for daytime
daytime_density <- nkde(
  bmr_roads_simple, 
  daytime_acc_projected,
  w = rep(1, nrow(daytime_acc_projected)),
  samples = daytime_acc_projected, 
  kernel_name = "quartic", 
  bw = 500, 
  method = "simple",
  verbose = FALSE
)

# Perform NKDE for nighttime
nighttime_density <- nkde(
  bmr_roads_simple, 
  nighttime_acc_projected,
  w = rep(1, nrow(nighttime_acc_projected)),
  samples = nighttime_acc_projected, 
  kernel_name = "quartic", 
  bw = 500, 
  method = "simple",
  verbose = FALSE
)

# Perform NKDE for weekdays
weekday_density <- nkde(
  bmr_roads_simple, 
  weekday_acc_projected,
  w = rep(1, nrow(weekday_acc_projected)),
  samples = weekday_acc_projected, 
  kernel_name = "quartic", 
  bw = 500, 
  method = "simple",
  verbose = FALSE
)

# Perform NKDE for weekends
weekend_density <- nkde(
  bmr_roads_simple, 
  weekend_acc_projected,
  w = rep(1, nrow(weekend_acc_projected)),
  samples = weekend_acc_projected, 
  kernel_name = "quartic", 
  bw = 500, 
  method = "simple",
  verbose = FALSE
)

Code

# Check if number of rows in density results matches the road geometries
if (length(daytime_density) != nrow(bmr_roads_simple)) {
  # Option 1: Pad the density with NA or 0 for roads without accidents
  daytime_density <- c(daytime_density, rep(NA, nrow(bmr_roads_simple) - length(daytime_density)))
}

# Repeat the same for nighttime, weekday, and weekend densities
if (length(nighttime_density) != nrow(bmr_roads_simple)) {
  nighttime_density <- c(nighttime_density, rep(NA, nrow(bmr_roads_simple) - length(nighttime_density)))
}

if (length(weekday_density) != nrow(bmr_roads_simple)) {
  weekday_density <- c(weekday_density, rep(NA, nrow(bmr_roads_simple) - length(weekday_density)))
}

if (length(weekend_density) != nrow(bmr_roads_simple)) {
  weekend_density <- c(weekend_density, rep(NA, nrow(bmr_roads_simple) - length(weekend_density)))
}

# Convert KDE results to sf objects
daytime_density_sf <- st_as_sf(data.frame(geometry = st_geometry(bmr_roads_simple), density = daytime_density))
nighttime_density_sf <- st_as_sf(data.frame(geometry = st_geometry(bmr_roads_simple), density = nighttime_density))
weekday_density_sf <- st_as_sf(data.frame(geometry = st_geometry(bmr_roads_simple), density = weekday_density))
weekend_density_sf <- st_as_sf(data.frame(geometry = st_geometry(bmr_roads_simple), density = weekend_density))

# Visualize KDE for daytime accidents
tm_shape(bmr_roads_simple) +
  tm_lines(col = "gray") +
  tm_shape(daytime_density_sf) +
  tm_lines(col = "density", palette = "Reds", lwd = 2, title.col = "Daytime Density") +
  tm_layout(title = "Accident Density during Daytime")

Code

# Visualize KDE for nighttime accidents
tm_shape(bmr_roads_simple) +
  tm_lines(col = "gray") +
  tm_shape(nighttime_density_sf) +
  tm_lines(col = "density", palette = "Blues", lwd = 2, title.col = "Nighttime Density") +
  tm_layout(title = "Accident Density during Nighttime")

Code

# Visualize KDE for weekday accidents
tm_shape(bmr_roads_simple) +
  tm_lines(col = "gray") +
  tm_shape(weekday_density_sf) +
  tm_lines(col = "density", palette = "Greens", lwd = 2, title.col = "Weekda Density") +
  tm_layout(title = "Accident Density during Weekdays")

Code

# Visualize KDE for weekend accidents
tm_shape(bmr_roads_simple) +
  tm_lines(col = "gray") +
  tm_shape(weekend_density_sf) +
  tm_lines(col = "density", palette = "Purples", lwd = 2, title.col = "Weekend Density") +
  tm_layout(title = "Accident Density during Weekends")

Code

tmap_mode("plot")